module compute
	use basicmodule
	use adjcvmod
	use dependence
	implicit none
	
	integer				:: nxsi

	real				:: tstconst
	real, allocatable	:: xsigrid(:)
	real, allocatable	:: taumsxgrid(:,:)

	
	contains
	
	
	subroutine Size_Dependence
		integer, parameter	:: ns=6, nsim=10000
		integer	:: l,j,t,ig
		real	:: stats(nsim),xsi,mu,sig,s
		real	:: msx(3)
		real	:: sizetab(nxsi,ns)
		call loadcvcoefs 
		sizetab=-1
		do ig=1,ns
			s=.1*(5.0/.1)**((ig-1.0)/(ns-1))
			do j=1,nxsi			
				xsi=0.03+(DGPxsimax-0.03)*(j-1)/(nxsi-1.0)
				msx=[0.0,1.0,xsi]
				call rnset(39184)
!				call depcorr(xsi,s*xsi);			
				do l=1,nsim
					call draw_dependence(msx,s*xsi)
					stats(l)=stat_adj(cvcoefs,gettest(msx))
				enddo
				sizetab(j,ig)=count(stats>1)/real(nsim)
				print *,sizetab(j,ig)
			enddo
		enddo
		print *,"sizetab"
		call mdisp(sizetab)
		print *,maxval(sizetab)
		call printtime
	end subroutine
	
	subroutine setnull(msx)
		real	:: msx(3)
		if(tstind==tstNyb) return
		selectcase(tstind)
			case(tstqnt)
				tstconst=msx(1)+(msx(2)/msx(3))*((-log(qGEV))**(-msx(3))-1)
			case(tstxsi)
				tstconst=msx(3)
			case(tstPar)
				msx(1)=msx(2)/msx(3)
			case(tstZpf)
				msx(3)=1
				msx(1)=msx(2)/msx(3)
		endselect
	end subroutine
	
	subroutine set_tststats()
		integer	:: l,j
		real	:: xsi
		real	:: msx(3),sig,mu,quants(nxsi)
		if(allocated(allstats)) deallocate(allstats,xsigrid)
		allocate(allstats(2,nsim,nxsi),xsigrid(nxsi))
	!call loadallstats;	
	!	do j=1,nxsi
	!		xsi=DGPxsimin()+(DGPxsimax-DGPxsimin())*(j-1)/(nxsi-1.0)
	!		if(abs(xsi)<1E-4) xsi=.001
	!		xsigrid(j)=xsi
	!	enddo
	!
	!return
		
		do j=1,nxsi
			if(taumodel) then
				mu=taumsxgrid(1,j)
				sig=taumsxgrid(2,j)
				xsi=taumsxgrid(3,j)
			else
				xsi=DGPxsimin()+(DGPxsimax-DGPxsimin())*(j-1)/(nxsi-1.0)
				mu=0;sig=1
			endif
			if(abs(xsi)<1E-4) xsi=.001
			msx=[mu,sig,xsi]
			call setnull(msx)
			xsigrid(j)=msx(3)
			do l=1,nsim
				call draw_Ys(msx,l)
				allstats(:,l,j)=gettest(msx)
			enddo
			quants(j)=count(allstats(1,:,j)>.7)/real(nsim)
		enddo
		call mdisp(quants)
		call printtime				
		call saveallstats
	end subroutine

	subroutine savexsicvs
		real	:: cvtab(nxsi)
		integer	:: j
		do j=1,nxsi
			cvtab(j)=quantile_s(allstats(1,:,j),0.95)
		enddo
		call savecsv(trim(dir)//"cvs_"//getfilesuffix()//".csv",xsigrid.cvr.cvtab)
	end subroutine
	
	subroutine saveallstats
		integer	:: myunit,ioerr
		open(file=trim(dir)//"allstats/stats"//getfilesuffix()//".fin",form="unformatted",access="sequential",newunit=myunit,iostat=ioerr)
		if(ioerr>0) then
			print *,"trouble opening file to save"
			stop
		endif
		write(myunit) allstats
		close(myunit)
	end subroutine		

	subroutine loadallstats
		integer	:: myunit,ioerr
		open(file=trim(dir)//"allstats/stats"//getfilesuffix()//".fin",form="unformatted",access="sequential",newunit=myunit,iostat=ioerr)
		if(ioerr>0) then
			print *,"trouble opening file to save"
			stop
		endif
		read(myunit) allstats
		close(myunit)
	end subroutine		
	
	
	function getllhess(msx0) result(hess)
		real	:: msx0(3),hess(3,3)
		integer	:: i,j
		real	:: msx(3),h(3,3),sumh(3,3)
		real, parameter	:: eps=epsilon(1.0)**.33
		sumh=0
		do j=1,n
			msx=msx0
			do i=1,3
				msx(i)=msx(i)+eps
				h(:,i)=getscore(msx,ys(:,j),yks(j))
				msx(i)=msx(i)-2*eps
				h(:,i)=(h(:,i)-getscore(msx,ys(:,j),yks(j)))/(2*eps)
				msx(i)=msx(i)+eps
			enddo
			sumh=sumh+h
		enddo
		hess=.5*(sumh+transpose(sumh))
	end function

	function gettest(msx0) result(val)
		use iercd_int
		real	:: val(2)
		real	:: msx0(3)
		integer	:: j
		real	:: msx(3),invhess(3,3)
		real	:: s(3),ss(3),lr
		msx=maxll(msx0,.false.)
		if(tstind.ne.tstNyb) then
			val(2)=msx(3)
			lr=getconstll(msx)
			lr=lr-getconstll(maxll(msx0,.true.))
			val=[lr,msx(3)]
			return
		endif
		invhess=-getllhess(msx)
		invhess=matinv(invhess)
		s=getdconstll(msx)
		if(sum(s*matmul(invhess,s))>1E-1) then
			if(msx(3)>-0.989) print *,"no interior maximum", msx
!			call mdisp(msx.cvr.s.cud.getllhess(msx))
			val=0
			return
		endif

		ss=0; val=0
		do j=1,n
			s=getscore(msx,Ys(:,j),yks(j))
			ss=ss+s
			val(1)=val(1)+sum(ss*matmul(invhess,ss))
		enddo
		val=[val(1)/n,msx(3)]
	end function

	subroutine obj(msx,iact,f,ierr)
		real	:: msx(3),f
		integer	:: iact
		logical	:: ierr
		ierr=.false.
		if(iact==0) then
			f=-getconstll(msx)
			if((f.ne.f) .or. f>1E10) ierr=.true.
			return
		endif
		if(tstind==tstNyb) return
		selectcase(tstind)
			case(tstxsi)
				f=msx(3)-tstconst
			case(tstqnt)
				f=msx(1)+(msx(2)/msx(3))*((-log(qGEV))**(-msx(3))-1)-tstconst
			case(tstPar)
				f=msx(1)*msx(3)-msx(2)
			case(tstZpf)
				if(iact==1) then
					f=msx(1)-msx(2)
				else
					f=msx(3)-1.0
				endif
		endselect
	end subroutine		

	subroutine dobj(msx,iact,g)
		real	:: msx(3),g(3)
		real	:: xsi
		integer	:: iact
		if(iact==0) then
			g=-getdconstll(msx)
			return
		endif
		if(tstind==tstNyb) return
		selectcase(tstind)
			case(tstxsi)
				g(1:2)=0
				g(3)=1
			case(tstqnt)
				g(1)=1
				g(2)=(1.0/msx(3))*((-log(qGEV))**(-msx(3))-1)
				xsi=msx(3)
				g(3)=(msx(2)*(-1 + (-Log(qGEV))**xsi - xsi*Log(-Log(qGEV))))/(xsi**2*(-Log(qGEV))**xsi)
			case(tstPar)
				g(1)=msx(3)
				g(2)=-1.0
				g(3)=msx(1)
			case(tstZpf)
				if(iact==1) then
					g=[1.0,-1.0,0.0]
				else
					g=[0.0,0.0,1.0]
				endif
		endselect
	end subroutine		

	function maxll(msx0,applyconst) result(msx)
		use nnlpg_int
		use cdgrd_int
		real		:: msx0(3),msx(3)
		logical		:: applyconst
		real		:: lbounds(3),ubounds(3)
		integer		:: m
		external	:: exobj, exdobj
		
		lbounds=[-1E3,1.000,-0.99]
		if(applyconst .and. tstind==tstPar) lbounds(3)=0.01
		ubounds=[1E3,1E3,10.0]
		if(applyconst) then
			m=merge(2,1,tstind==tstzpf)
		else
			m=0
		endif
		call erset(0,0,0)
		call nnlpg(exobj,exdobj,m,m,0,lbounds,ubounds,msx,xguess=msx0+[0,1,0],iprint=0,maxitn=500)
		msx=msx-[0,1,0]
	end function
	

end module


	subroutine exobj(msx,iact,f,ierr)
		use compute
		real	:: msx(3),f
		integer	:: iact
		logical	:: ierr	
		call obj(msx-[0,1,0],iact,f,ierr)
	end subroutine
	
	subroutine exdobj(msx,iact,g)
		use compute
		real	:: msx(3),g(3)
		integer	:: iact
		call dobj(msx-[0,1,0],iact,g)
	end subroutine
	